Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 26 January 2013 (MN WT^ style file v2.2) 



Three Body Resonance Overlap in Closely Spaced Multiple 
Planet Systems 

Alice C. Quillen 

Department of Physics and Astronomy, University of Rochester, Rochester, NY 14627, USA 



26 January 2013 



o 

(N 



Oh! 

O 



(N 
> 

in 

o 

\6 
o 



ABSTRACT 

We compute the strengths of zero-th order (in eccentricity) three-body resonances 
for a co-planar and low eccentricity multiple planet system. In a numerical integra- 
tion we illustrate that slowly moving Laplace angles are matched by variations in 
semi-major axes among three bodies with the outer two bodies moving in the same 
direction and the inner one moving in the opposite direction, as would be expected 
from the two quantities that are conserved in the three-body resonance. A resonance 
overlap criterion is derived for the closely and uniformly spaced, equal mass system 
with three-body resonances overlapping when interplanetary separation is less than 
an order unity factor times the planet mass to the one quarter power. We find that 
three-body resonances are sufficiently dense to account for wander in semi-major axis 
seen in numerical integrations of closely spaced systems and they are likely the cause of 
instability of these systems. For interplanetary separations outside the overlap region, 
stability timescales significantly increase. Crudely estimated diffusion coefficients in 
eccentricity and semi-major axis depend on a high power of planet mass and inter- 
planetary spacing. An exponential dependence previously fit to stability or crossing 
timescales is likely due to the limited range of parameters and times possible in inte- 
gration and the strong power law dependence of the diffusion rates on these quantities. 



INTRODUCTION 



The stability of multiple planet systems has been long 
been a matter of interest as it concerns the long term 
stability of the Solar system. Poincare noticed that per- 
turbation techniques involved singularities or small divi- 
sors that prevented solution by convergent series. Recent 
numerical explorations suggest that the giant planets in 
our solar system were originally located in a more com- 
pact location and experienced a subsequent planet-planet 
scattering event (w ithin the context of the "Nice model"; 
iTsiganis et al.|[2005h . All extrasolar planetary s ystems may 
exper i ence epochs of dynamically instability (e.g. Ford et al.' 
20011: iBarnes fc G rccnbcre 2006; Chatterieo ct al. 20^ 



ies ar e sufficiently distant f r om ea c h oth er | Chambers et alj 
19961: iDuncan fc Lissaueil Il997l: |M )er & Quillen 2001 



Chatteriee et al.ll2008l : ISmith fc Li ssaucr 2009). For an inte- 



gration begun with all bodies at zero eccentricity and incli- 
nation, the integration time until one body crosses the orbit 
of an other body is known as a crossing timescale, tc- We 
refer to the idealized problem studied by Chamber s et al.l 
(119961) with equal mass planets and with the semi-major 
axis of each consecutive planet set from that of the previous 
one using a constant 5, 



an+i = (1 + 5)a„ 



(1) 



Thomnies et al.l 



.2 008: Go zdziewski fc Migaszewskil 

20091: I Raymond et al.i i2009al: iKopparapu fc Barnes! 



Fabrvckv fc Murrav-Clay||2010l ) 



2008 
201C 



Numerical integrations show that a system of two 
planets on initially zero-inclination and eccentricity or- 
bits about a star never experience mutual close en- 
counters if the initial scmimajor axis separation i s suf- 
ficiently large jMarchal fc Bozis 1982; Gladmaj 1 19931 : 
iBarnes fc Greenberd |2007| ) . (also see iMardlinel 12008 ). Sta- 
bility of multiple planet systems is often discussed in 
terms of this limit which has been called "Hill stabil- 
ity" (e.g. , 'Marchal fc Bozis" 198^; IBarnes fc Gree nberg"2007|: 
[Raymond ct al. 2009b). Systems with multiple planets or 
satellites are stable over long periods of time if the bod- 



The planets have mass ratio with respect to the central star 
m = rrip/Mf When the planet number is greater than five, 
the c rossing timescale is ins ensitive to the number of bodies 
(e.g., IChanibers et al.lll996l ). 

A number of stu dies have fit power laws to the crossing 
or stability timescale. iDuncan fc Lissaueil l| 19971 ) found that 
the crossing time is sensitive to the mass of the satellites 
wit h logtc ~ otm -\- P where a is a slope and /3 is an off- 
set, [chambers et al. I Il996) foun d logic ~ a(5m-^'"') + /3, 



whereas Smith fc Lissauerl lj2009l ) fit log tc 



t(5f 



-1/3n 



Here the parameters a, /3 are fit to the numerically measured 
crossing timescales and are not identical in each setting. 
These numerical studies integrated for timescales between 
100 and 10^ orbital periods of the innermost orbiting body. 
The interplanetary separation ranged from 1 to about 10 
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mutual Hill radii and the mass ratio ranged from m ~ 10~ 
to 10~^. The trends in the numerically measured crossing 
timescales are currently lacking an explanation. 

The power law forms for the crossing timescales 
can be re-written tc oc exp(m"5^) for exponents a, 7 
(IFaber fc Quilled [2007! '). The exponential form is reminis- 
cent of the Nekhoroshev theorem ( Nekhoroshe v 



of Arnold diffusion (| Arnold! 1 196'3 ) in the context of weakly 
perturbed Hamiltonian systems. There are subtle and deep 
connections between the Nekhoroshev theor em and Arnold 
ditfusion (e.g, as explored and discus sed by IChirikov|[l979l : 
iLochak fc Neishtadt) [l993 ; [Lochak 199^ 1. Arnold diffusion 
takes place on exponentially long timescales whereas sys- 
tems with resonances that are fully overlapped can have rel- 
atively faster diffusion rates. This led MorbidcUi & Frocschlc 
(|l995l ) to suggest that systems with sparse and weak 
resonances might exhibit diffusion on exponentially long 
timescales whereas those affected by multiple and overlap- 
ping resonances would diffuse at a rate that d epends on 
a power law of time (also see IGuzzo et al.ll2002l ). A num- 
ber of studies have numerically measured a power law rela- 
tion between Lyapunov and instability or cro s sing timescales 
jLecar et al.l 1993: Levison fc Duncad 19931: Murison et al.| 



I994I : iMorbidcUi fc Froeschle 'l995': Mikkola & Tanikawa! 
2OO7I : lUrmin skv fc Hcggc 2009; Shcvchcnko 2010) suggest- 
ing that the stability or crossing timescales are driven 
by chaotic diffusion or Hamiltonian intermittency. Because 
there is a power law relation between Lyapunov and crossing 
timescales, we would expect that the diffusion rate would be 
a power law function of the perturbation parameters (such 
as planet mass and planetary separation) rather than an ex- 
ponential function of them as commonly fit. This apparent 
contradiction has been discussed in terms of two regimes for 
the dynamics, a weakly perturbed "Nekhoroshev" regime 
exhibiting diffusion over exponentially long timescales and 
a resonant overlap power law regime exhibiting diffusion 
at a rate that de pends on a power of the perturbations 
(|Morbidelli fc Froe schle 1995,). 

In this study we consider the role of three-body reso- 
nance s in the idealized setting studied by [Chambers et al.l 
(|l996t ) of a closely and uniformly spaced equal mass co- 
planar system initially in nearly circular orbits. Our goal 
is to understand the dynamics of the idealized coplanar, 
closely and uniformly spaced equal mass multiple planet 
system sufficiently well that we can identify the source of 
instabilities in multiple planet systems. Previous studies 
of three-body resonances have primarily focused on set- 
tings where one of the bodies is small, suc h as an aster- 
oid perturbed by both Jupiter and Saturn ([Murrav et al.l 
1 1998 : iNesvornv fc Morbidelli|[l999l : iGuzzolbOOsf T but also in- 
clude the early study of the Laplace resonance by lAksned 
(| 19881 ). While three-body resonances are weak there are 
more of them than two-body resonances so they can be a 
source of chaotic beha vior causing slow diffusion in eccen- 
tricity and inclina,tion (INesvornv fc Morbidellil 1 19981 . 1 19991 : 
IGuzzo et al.|[2003 : lGuzzdl20'05l )7 Three body resonances may 
be important in extra solar multiple planet systems. The 
long-term stability of extra solar multiple planet systems 
may be influence d by a net of low order two and three- 
body resonances (^Gozdziewski fc Migasze^^ |2008| . l2009l : 
iFabrv ckv fc Murray-Clay 2010). 

We first estimate the strength and libration frequency in 



zero-th order (in eccentricity) three-body resonances. Con- 
served quantities for them are also estimated so that their 
signature in numerical integrations can be identified. Using 
estimated numbers and widths of the three-body resonances 
we derive a resonance overlap criterion. In the final section 
we discuss the role of three-body resonances in causing in- 
stability in multiple planet systems and whether they may 
eventually provide an explanation for the exponential forms 
fit to their crossing timescales. 



2 HAMILTONIAN FOR A MULTIPLE 
PLANET SYSTEM 

2.1 Non-interacting System 

The Hamiltonian for A'^ non-interacting massive bodies or- 
biting a star (and so feeling gravity only from the star) can 
be written 



Ho = 



JV 

E 



"2Af 



(2) 



where rrij is the mass of the j-th body (or planet) divided 
by the mass of the star, M, . Here we have ignored the mo- 
tion of the star and have put the above Hamiltonian in 
units such that GM« = 1 where G is the gravitational con- 
stant. Here the Poincare coordinate A, 



the semi-major axis of the j-th planet is aj 
coordinate is conjugate to the mean longitude, 
j-th body. The mean longitude Xj = Mj + zuj where Mj 
is the mean anomaly and zui is the longitude of pericenter 
of the j-th body. We may also use the Poincare coordinate 

r, = 



mjy^/Oj, where 
. This Poincare 
A, of the 



/5J(1 - yi^) 



Ltjej/2 where ej is the 



j-th body's eccentricity. This coordinate is conjugate to the 
angle 7j = —'ujj. We will restrict our system so that all plan- 
ets are orbiting in the same plane and so will begin by ig- 
noring the Poincare coordinates associated with inclination 
and the longitude of the ascending node. Each Poincare mo- 
menta contains a factor of a planet's mass. Some studies of 
three-body resonances have focused on the problem of a low 
mass object in the presence of two p lanets (e.g., an asteroid 
perturbed by Jupiter an d Saturn; INesvornv fc Morbidellil 
1 19981 : iMurrav et al.l Il998l ) and so have removed the mass 
from the Poincare momenta associated with the low mass 
object. 



2.2 Interactions 

We now consider the gravitational interactions between the 
planets. Here we consider only the direct term and ignore 
the indirect terms. The Hamiltonian can be written 



(3) 



where the interaction term, Hint, is a sum of direct interac- 
tion terms. Hint = ^^-^ . (r;, rj), and 



If the planets are in nearly circular orbits 



(4) 



(5) 



3jV^ + "ij ~ 2a,i cos(A, - Xj) 
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Figure 1. Approximations to tlie Laplace coefficient. Plotted as 
points are the Laplace coefficient h^^^^{a) as a function of the in- 
teger q for four separations 5 = a^^ — 1 = 0.01, 0.05, 0.1, and 0.2. 
Overplotted as thick brown lines for each 5 value is the function 
0.5| ln(<5)| exp(— g(5) and as thin black lines 0.5j ln(<5)| exp[— g<5(l — 
0.002 In 5)]. 



and we have assumed that aj > at and 



a 



13 



We can expand in Fourier components 



q=0 

with coefficients 



i{q\i — qXj) 



where b')^j^{a) is a Laplace coefficient. 



cos{q(j))d(j) 



(1 + q2 



2q cos ( 



(6) 



(7) 



(8) 



(9) 



Laplace coefficients are the Fourier coefficients of twice the 
function f{4>) = (l + a^ — 2q: cos 0)"^. As this function is lo- 
cally analytic the Fourier coefficients decay rapidly at large 
q and the rate of decay is related to the width of analyt- 
ical continuation in the complex plane. This function can 
be analytically continued on the complex plane in the re- 
gion a < |z| < with f{z) = (1 + - a{z + z~^)~" = 
i 6"(Qf)2:". The Cauchy root test for convergence 

implies that in the limit of large n that |6"(q:)| < a" and 
so the Fourier coefficients decay rapidly. We approximate 
Laplace coefficient with the function 



(10) 



b[)[{a) ~ 0.5lln5|exp(-<51q|) 

where 5 is the interplanetary separation (equation [T]) and 
5 = — 1 ~ 1 — a. This Laplace coefficient diverges log- 
arithmically for small S (or for a near 1) and drops expo- 
nentially at large q. In Figure [T] we graphically show this 
approximation for the Laplace coefficient for S in the range 
0.2 to 0.01. 

We can write the interaction term as a function of 
Poincare coordinates. 



^1/2 



A2m2 



(11) 



Indirect terms can be neglected here because they 
only contribute a single zero-th order Fourier compo- 
nent, that with q = 1 (e.g., see appendix Table B.2 by 
iMurrav fc Dermott|[T999l '). 



3 THREE BODY RESONANCES 

We consider the possibility that the system is not in any 
two body resonances where pru ~ qrij with integers p, q, 
but might be in a three body resonance. Here Ui is the mean 
motion of the i-th body. A Laplace relation exists between 
three orbiting bodies if the frequency 



pui — (p + q)nj + quk ~ 0. 



(12) 



Integrating the previous equation (and assuming that the 
precession rates are slow) we find that the angle 



= pXi — {p + q)\j + q\k ^ constant. 



(13) 



This angle librates about a particular value (often or vr) 
when in a three-body resonance. The period, T, corresponds 
to the time between successive repetitions of the initial con- 
figuration 



T 
2^ 



p + q 

rii — Uk 



(14) 



We try to maintain the ordering < < and > 
rij > rife. 

We explore how 3-body interaction terms involving an- 
gles such as given in equation p3p can be constructed from 
individual 2-body interaction terms. A similar procedure has 
been used be f ore to estimate 3-body resonance strengths 
dAksnej [l988l: iNesvornv fc Morbidellii [l998l: iMurrav et all 
19981: iNesvornv fc MorbideUil Il999l : IGuzzoI |2005| ) (also see 
Chirikovlll979i r Our procedure is to carry out a first order 
canonical transformation that is designed to remove the two 
2-body interaction terms. The procedure is d escribed, for 
example, in section 4.1 bv lFerraz-Meiiol l|200'i1) . 



3.1 Canonical transformation removing first order 
(in mass) two-body interaction terms 

We begin with the Hamiltonian for three bodies, lacking in- 
direct terms and with two Fourier components from two sep- 
arate two-body interaction terms. We first consider zero-th 
order components (in eccentricity) so we need only consider 
the Poincare coordinates A, A where the vectors refer to the 
the coordinates and momenta for three planets. We chose 
two components with arguments (angles) whose difference 
is equal to a three-body Laplace angle (equation [13} 

p{\i — Xj) — q{Xj — Afe) = pXi - {p + q)Xj + qXk. 

Our Hamiltonian with these two components 



l = i,j,k 



2Af 



(15) 

{pXi - pXj) + Wjk,q cos{qXj - qXk) 
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with the functions 



A? 



"1/2 



(16) 



as discussed in the previous section. We have used the no- 
tation 



Oljk 



nijAi 



It is useful to compute some derivatives 



2mim 
"T^Af 



dAi 



[l + a^jDc]bl.^{aij) 



dW^ 



dAi 



2m, m? 
— -I 



rfAfe 
with Da 



[l + Qjfe-D„]bJ/2(ajfc) 
We can use the shorthand 



(17) 



(18) 



and similarly with indices jk. 

We use a generating function with new momenta A' and 
old coordinates A 



F2{A',X) 



E a;a. - 



l—i,j,k 

W'l. 



wi 

— ^ sin(pAij) 
pn-. 



qn 



sin(gA 



jk) 



(19) 



'jk 



to generate a canonical transformation. Here is a func- 

tion of A^ and A^ and similarly n'^j . This canonical transfor- 
mation is designed to remove the two perturbation terms in 
the Hamiltonian to first order in the planet masses. Deriva- 
tives of the generating function give us new coordinates in 
terms of the old ones 



Ai 



Afc 



A' 



W 



3(pAij) 



(20) 



W 

a; + ^ cos(pA, 



W'l. 

" jk.q 



cos(gAjfc) 



\l I ^' jk,q 

J^k H ; — 



cos{qXjk) 



X[ = X^ + 



A' = 



dW 



dn', p 

aA; png 9A^ pn[. 

dAr pn[] 



^k 



+ 



Xk 



9w;k,q 



pn., 
1 



9^ 



dK mfk 

It is useful to relate 



dn'k W',k,q ^ dW'.k.q 



sin(pAij) 
sin(pAij) 
sin(gAjfc) 
sin(ijAjfc) 



Ml 



3n4 



(21) 



and similarly for the other bodies. 

We replace our old coordinates and momenta with new 
ones in the Hamiltonian finding that in the new coordinates 
the first order two body terms have been removed by the 
transformation. We expand to second order in the masses 
and find that the Hamiltonian has gained second order 
terms. During this procedure we drop terms that depend on 
cos^(pAij), or sin^(pAij) and cos^(gAjfc), or sin^(gAjfe) while 
keeping those containing the products cos(pAij) cos(gAifc) 
and sin(pAij) sin(qAjfc). We rewrite these products in terms 
of the sum and difference of the angles and discard the term 
that contains the sum of the angles so as to retain only the 
term that depends on the Laplace angle qXi — {p+q)Xj +qXk- 
When the Laplace angle is slowly varying (and the system 
is near a three-body resonance) the other terms are rapidly 
varying and so can be neglected. 



3(pA:-(p + g)A;+gA',)(22) 



with 













1 


2 \2nr.n',^ 








+ (4+ 


pn'j 


I'^'jk 



qn'. 



+ 



P^?k 



&^/2(a-. 



bl/2(^'^j)°''jkDabl^^{a'jt 



(oi'jk) 
(23) 



Our procedure using a canonical transformation should give 
a resona nce term equivalent to the zero-th order term de- 
rived bv I Aksnesi ()1988i) using Lagrange's equations, though 
it is not easy to check because of the differences in notation. 
The resonance term is zero-th order in planet eccentricity 
but second order in the planet masses. Hereafter we drop 
the primes in the coordinates. 

The full Hamiltonian contains additional two-body 
terms however the angles involved are expected to vary 
quickly compared to (f>. Neglecting fast angles is equivalent 
to averaging over them. Equivalently as long as there are 
no combinations that yield slow angles, the other interac- 
tion terms may be removed using near identity canonical 
transformations similar to that used above. 



3.2 Width, libration frequency and conserved 
quantities for the zero-th order three-body 
resonance 

The Hamiltonian (equation I22|l that contains a three-body 
term can be used to estimate the width and timescales in a 
Laplace resonance. We can perform a canonical transforma- 
tion to reduce the dimension of the problem. Consider the 
generating function 

F2{X, J) = (pA, - (p + q)Xj + qXk)J + XjJj + XkJk (24) 
leading to new angles (</>, X'j, Aj.) 
(f) = {pXt - {p + q)Xj + qXk) 
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and new momenta (J, Jj, Jk) such that 

pj = Ai 
-(p + q)J + Jj = Aj 
qj + Jk = Afc . 



(25) 



(26) 



After the transformation, the new Hamiltonian (using 
equation 1221) is 



B = z • no (34) 

setting distance to resonance, where the vector of integers 

z={p,-{p + q),q). (35) 

On resonance B ~ 0. The coefficient A depends approx- 
imately on the magnitude of the vector z with \A\ ~ 
lzlV(m-a«)- 

With a shift in the momentum. Is — I + B/A, the 
Hamiltonian (equation I32p can be written so as to remove 
the term that is proportional to Is, 



K{J,Jj,Jk;<p,\j,Xk) = 



2p2J2 2(Jj - (p + g)J)2 

+ epq COS (p. (27) 



ml 



2{Jk + qjy 



The new Hamiltonian only depends on the angle and does 
not depend on the two longitudes \j, Xk (that are unchanged 
by the canonical transformation) so the two conjugate mo- 
menta Jj , Jk are conserved quantities. Our conserved quan- 
tities can be written as 



pjj = pAj + {p + q)Ai 
pJk = pAk-qAi. 



(28) 



Differentiating these conserved quantities with respect to 
the semi-major axes we find that small changes 

dak _ m^q f 2±\ 
dai ruk p \ai J 

daj _ rrii (p + g) ( aj\^^^ 
dai Ti 



P 



(29) 



These derivatives imply that a small change in semi-major 
axis by one body will be mirrored by changes in semi-major 
axis of the two other bodies. The signs imply that the outer 
two bodies move in the same direction but the middle one 
moves in the opposite direction. The middle body is ex- 
pected to move more than the outer two as we expect p + q 
is greater than p and q. This behavior can be seen in particle 
integrations as we will discuss below. 

We can expand the momentum J about an initial value. 
Consider initial values for Aio,Aj,o,Afco corresponding to 
initial value Jo, conserved quantities Jjo,Jko, initial semi- 
major axes and mean motions 



ao = (flio, fljo, afco) no = (nio,njo,nko)- 
We define 
J = Jo+I, 



(30) 



(31) 



and expand the Hamiltonian (equation I27|) about Jq. To 
second order in I 



K((j}, I) = + BI + epq cos (j) + constant 



(32) 



where the constant contains terms that depend on our con- 
served quantities {Jj,Jk) and Jo. The coefficients 



B = puio - {p + q)njo + quko 

„2 / I \2 2 



A = 



(p + q) ^ r 



(33) 



The coefficient epq is also evaluated at ao. We can think of 
the coefficient B as a product 



K{(f>, Is) = A-^ + epq cos (j) + constant. 



(36) 



We estimate the width of the resonance in momentum is 

(37) 



A/~2a/^ 



corresponding to a resonant width in Poincaro momentum 
Ai of 



AA. 



(38) 



(using equation [26} and a width in terms of semi-major axis 
of the innermost body 



4p 2epqai 



A 



(39) 



A jump across resonance would give a change Aa^ for the 
inner body. Using equations (|29|l changes in semi-major axis 
for the other two bodies can be estimated from that of the 
inner one. 

The libration frequency in the resonance is 



(40) 



For a system initially with J = to be in resonance we 
require that the shift, B/A, (relating I and Is) is smaller 
than half the resonance width (A//2 in eauation l37|l . Using 
this condition and equations (|37p and (|40p we require jBj < 
\piupq to be near or in resonance or using equation (|34p 



z ■ no 



< V2 < 



(41) 



3.3 Estimates of three-body resonance strengths 
and frequencies for closely and evenly spaced 
equal mass multiple planet systems 

We consider the strength of the various terms contributing 
to the Laplace resonance strength, epq in equation for a 
closely and evenly spaced equal mass system. For the equally 
spaced system (represented by equation[l} 5ij ~ &jk - In this 
setting the only Laplace angles that can be nearly fixed have 
p about the same size as q. We work in units of the mean mo- 
tion and semi-major axis of the innermost body involved in 
the three-body resonance. Differences in the mean motions 
can be approximated as 



where 



(42) 



(43) 
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As the Laplace coefficient b^-^J^ (a) can be approximated given 
in equation (|10|l . the derivatives of the Laplace coefficients 
can be approximated as 

Dc.h'^^'j^ia) ~ 0.5(5"' + \q\n5\) cM-5\q\) (44) 

(using Da = —js)- Using these approximations and assum- 
ing equal masses, we find that the interaction term strength 
in equation (|23|l is approximately 



€pg ~ [S^'^ilnSf + 0.5S'^ {2 + {p + q)5\\n5)\)\\nS\ 



X exp(— (5(p + q)). 



(45) 



where we have assumed integers p, q are positive. The terms 
all have the same sign and in most cases the first term dom- 
inates so we can approximate the interaction strength as 



m S (In 5) exp{—S{p + q)). 



(46) 



Note that epq is positive, so we would expect libration 
around Laplace angle of zero or 



= pXi - {p + q)\j + q\k ~ 0. 



(47) 



For p ~ q (corresponding to being near resonance for 
small S) we estimate that |j4j ~ 20p^/m (using equation 
I33p and from equation (|40p and equation (|46p the libration 
frequency in the resonance is approximately 



• 4mpS ""^1 ln<5| exp(— 5p). 



(48) 



The width of the resonance (in terms of momentum Ai ; 
equations 1381 [5^ gives a width in semi- major axis similar in 
size or 



Aai ~ mS ^\\nS\ 



(49) 



where we have assumed pS < 1 and so neglected the expo- 
nential. The above equation should give the size of jumps 
across resonance for the inner body. 

Likely equation (|49|l somewhat underestimates the res- 
onance widths as we have not taken into account all terms 
in equation (|23|l . We have also neglected the dependence 



of epq on the momentum J in our estimates of copq (equa- 
tion 30} and Aai (equation I39p . The conserved quantities in 
resonance imply that 5ij increases when Sjk decreases and 
vice-versa, so epq may not be strongly dependent on J. 

Because the Laplace coefficients drop exponentially 
with separation 5 (see Figure [1} or the distance between the 
planets, resonances between the first, second and fourth bod- 
ies or other non-consecutive combinations should be much 
weaker than those involving three consecutive bodies. How- 
ever if the masses of the bodies difler then three-body res- 
onances involving non-consecutive triplets could be impor- 
tant. 



3.3.1 First order resonances 

Zero-th order (in eccentricity) resonances do not influence 
the Poincare variable associated with eccentricity so they 
should not affect planet eccentricities. Fi rst order three -body 
resonances (as previously considered bv lAksne3ll988l ) have 
interaction terms in the form 



V^y^ cos(pAi - {p + q- l)Aj + qXk - zui) 



(50) 



with a single longitude of pericenter zui. The longitude of 
pericenter can be that of any of the three planets involved 



in the resonance, so that I can be equivalent to i,j or k. Be- 
cause the interaction term contains a longitude of pericenter 
it can affect a planet's eccentricity. The strength rj can be 
estimated in the same way as we have estimated the zero-th 
order three-body resonance strengths but instead of carry- 
ing out a transformation with two zero-th order two-body 
interaction Fourier terms, one begins with a zero-th order 
and a first order two-body Fourier term. 

When expanded to first order in planet eccentricity and 
inclination the two-body interaction terms (equation^)) gain 
Fourier components (that would be added to equation [7| 

oo 

[Vrj,q cos(gAj + (1 - q)Xi - ^0 + 

q— — oo 

Vij,q cos(gAj + (1 - q)\i - -c^j)] (51) 



where 



e^ f 27 (aij,q) ^ ^ (, Afj /27(Qy, 

J A 



^ij,q = SJ^ej /31 (Uij ,q)^ jj^ /31 (Uij , q) 



and coefficients 



/27(a,g) = -[-2q-aDa\b[%{a) 



/3i(a,q) = -hl + 2g + QD, 



-'1/2 



(a) 



(52) 



(53) 



(equation 6.107 lMurrav fc Dermottlll999l : also see Tables B.4 
and B.7). An approximation to these coefficients is — /27 ~ 
/31 ~ 0.5(5"' + jgln5|)exp(-5|g|) shown in Figure^ Also 

First order and zero-th order terms when combined 
form a term with a three body argument similar to that 
shown in equation (|50p . We list below on the left the two 
Fourier components that when combined give the argument 
on the right; 



y^^j,-p{o^^3Wjk,q{Ol■jk) 

Vij- -p(a»i)Wjfc,,(ajfc) 
Vjk,q{ctjk)Wij,p{aij) 

Vj'k,q{ajk)Wij,p{aij) 



[p + - (p + q)\j + q\k - T^i 
{p + l)Xi — (p + q)\j + qXk — 



p^i 
pK 



(p + q - l)Aj + q\k 
{P + 



+ qXk — u^k. 

(54) 



The longitude of pericenter in the argument determines the 
eccentricity related Poincare coordinate. For example, if the 
argument contains vji then the three-body term contains a 
factor of and similarly for bodies j, k. 

Whereas the zero-th order three-body resonances in- 
volve two W coefficients, the first order three-body reso- 
nances involve a single W and a single V coefficient. The Vij 
coefficients are approximately derivatives of the Wij coeffi- 
cients. The strength, r;, of the first order three-body term 
we expect is larger than epq by a factor of because it 
would involve an extra derivative of the Laplace coefficient. 
However the dependence on the Poincare coordinate F asso- 
ciated with the eccentricity can reduce the strength of the 
resonance. The first few conjunctions of an initially zero ec- 
centricity system lead to planet eccentricities of order a few 
times mS~^ (e.g., equation 10.57 Murray & Dermott 1999). 
This suggests that the ratio of the first to zero-th order 
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Figure 2. Approximations to the coefficients /27,/3i- Piotted as 
points are tiie coefficients — /27(a, ?) and /3i(a,g) (equations I53I I 
as a function of tiie integer q for four separations 5 = — 1 = 
0.01,0.05,0.1, and 0.2. Overpiotted as thin black lines for each S 
value is the function 0.5(5"^ + |gln<5|) exp(— These functions 
are approximately equivalent to Dctb^^,-'^{a). 



resonance strengths is approximately rjVj/epq ~ m5~^ . It 
may be convenient to define 5h = S/vh where cc m}^^ 
is the Hill radius of the planet. The ratio of resonance 
strengths is then approximately 35^'', or the interplanetary 
separation in units of the planet's Hill radius. For e ccentric- 
ities above a critical value (e.g., set dimensionally; iQuillenI 
l2006f ) the resonant width and libration frequency depends 
on the square root of r]Vj. As &h is in t he range ^ 2-10 
for the systems studied num erically (e.g., [Chambers et al.l 
1 19961 : ISmith fc Lissaueill2009l ') we expect that first order res- 
onances are initially a few times weaker than the zero-th 
order ones. However, a system evolving in multiple three- 
body zero-th order resonances can cross first order reso- 
nances causing variations in planet eccentricities. Because 
of the number of possible angular combinations there are 
more first order resonances than zero-th order ones. 



4 THREE-BODY RESONANCES AS SEEN IN 
A NUMERICAL INTEGRATION 

In Figure[3]we show a numerical integration of 5 equal mass 
bodies initially in a coplanar circular orbits about a central 
star. The ratio of the planet masses to that of the central 
star is m = 10~^. The initial separation between the bodies 
is given by 5 = 0.11247 using equation ((T|. The numerical 
integration was done using the hybrid algorithm of the code 
Mercury version 6.2 (Chambers 1999). Time is given in units 
of the initial rotation period of the inner body. Distances are 
given in units of the inner body's initial semi-major axis. 
This numerical integration was chosen to illustrate phenom- 
ena associated with three-body resonances and we will use it 
to check predicted sizescales for them. For these parameters 
^m"^/^ = 5.22 and 5m~^''* = 2.0. In te rms of the mutual 
Hill r adius (as defined by equation 1 by ISmith fc Lissaueil 
bOOgl ) (5 = 5.6r,ng, placing it in the middle of the regime 
studied by ISmith fc Lissaueil (|2009l ) though they primar- 
ily considered planets lower in mass by a factor of 3. This 
separation places our integration at larger separations than 



the mean explored by [Chambers et alj l|l996t) (on the right 
hand side of the top panel of their Figure 3) with a crossing 
timescale of about 10^ orbital periods. 

In Figure[3^ we show the semi-major axes of the 5 bod- 
ies. Variations in semi-major axis often involve similar mo- 
tions for three consecutive bodies, with the inner and outer 
ones (of this triplet) moving together and the central one of 
the triplet moving in the opposite direction. This is expected 
as in a three-body resonance there are two conserved quan- 
tities (equations I28p that relate variations in semi-major 
axes between the three bodies (equations [29]). Figure [3}3 
shows the Laplace angles = p\i — (p -I- q)\j + qXk for 
p = 5,p + q — ll,q — 6. We can write (j> = z ■ \ with 
z = (5,-11,6). The angle is plotted for the inner three 
consecutive bodies (bottom panel), the middle three con- 
secutive bodies (middle panel) and the outer three consecu- 
tive bodies (top panel), respectively. Figure[3J: is similar but 
for the resonances with p = 6,q — 7. We see from Figures 
^jp,c that there are times when specific Laplace angles vary 
slowly or vibrate about or tt. The system is affected by 
more than one Laplace resonance. At some times the inner 
three bodies move together, and at other times the middle 
or outer three move together. For example, the variations in 
the outer three planets at t = 19000 periods are likely due 
to the z = (6,-13,7) resonance involving the outer three 
bodies. Variations in the inner three planets at t = 9000 
are likely due to the z = (5,-11,6) resonance involving 
the inner three bodies. When the Laplace angle ceases to 
circulate and librates about or tt we see that variations 
in the three planets involved in the Laplace resonance are 
related, with the outer two increasing or decreasing in semi- 
major and the middle one moving in the opposite direction. 
The middle body experiences larger variations in semi-major 
axis as would be expected from the conserved quantities (see 
equations 1281 and I29|l . We have checked that the conserved 
quantities in these resonances do not make large variations 
when the resonant angle is not rapidly circulating. However, 
conserved quantities associated with one resonance can vary 
while a different resonance is affecting the system. 

For m = 10~^ and 5 — 0.11247 a jump across resonance 
should give a change in semi-major axis for the inner body 
(using equation I49p of approximately 2 x 10~*. Jumps seen 
in the simulation (Figure[3^) are a factor of a few larger than 
this. We can consider this moderately reasonable agreement 
as we have only made rough estimates of the resonance prop- 
erties. The libration frequency in the resonance (computing 
ujpq using equation I48p is approximately ujpq ~ 10"^ corre- 
sponding to a period of T = -22- ~ 600 years. This period 
is short enough that the slow variations in angle in Figures 
[3^,b can be attributed to the three-body resonances. 

The simulation shown in Figure [3] is affected simultane- 
ously by more than one zero-th order three-body resonance 
(here z = (6, —13, 7) and (5, —11, 6) resonances and for three 
possible consecutive triplets of planets) suggesting that they 
are dense and wide enough that the three-body resonances 
overlap. If so then random variations in semi-major axis can 
be attributed to chaotic behavior associated with multiple 
resonances. 

For the same numerical integration we also show the 
eccentricity evolution in Figure [4^. This figure also plots 
angles = z • A for z = (9,-14,4) and z = (2,-10,9). 
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These angles are not in the form z = {p,—(p + q), q) as the 
sum of the indices in the vector z is not zero. First order (in 
eccentricity) thee-body resonances involve a single planet's 
longitude of pericenter, for example, the angle could be one 
of the following 



(p + l)Ai — (p + q)\j + q\k — T^i 
p\i — (p + g — l)Aj + q\k — vjj 

p\i — (p + + {q+ l)Afc — -ajfe. 



(55) 



These angles arise from combining a zero-th order two-body 
term with a first order two-body term. As the precession 
rates, ro, are slow compared to the mean motions, we have 
plotted angles omitting a longitude of pericenter. We have 
examined similar plots containing all of the above possible 
angular combinations and found that they are similar to but 
noisier than Figures |4}3,c. Since the planet eccentricities are 
low, small variations in the orbits can give large changes in 
the computed longitudes of pericenter. 

From Figure|4j3,c we see that our numerical integrations 
also exhibit slow angles associated with first order three- 
body resonances. When first order resonances are crossed 
we expect small changes in planet eccentricity. For exam- 
ple at t ~ 8000 the inner three planets are affected by the 
z — (9,-14,4) resonance leading to an increase in eccen- 
tricity in these three planets. At t ~ 10, 000 years the outer 
three planets are affected by the z = (2, —10,9) resonance 
leading to an increase in the eccentricity of the fourth planet. 
Jumps in eccentricity seem similar in size to those of jumps 
in semi-major axis though we expected them to be a few 
times smaller based on the discussion in section 3.3.1. A 
comparison between the angles shown in Figure [Sja and Fig- 
ure |3J3 shows that they are very similar; likewise for Figure 
and Figure [3J;. Hence the zero-th order resonances (with 
angles shown in Figure [3l3,c) overlap the first order reso- 
nances (with angles shown in Figure [3)d,c). 

Once the eccentricity of a planet is increased secular 
perturbations cause oscillations in the eccentricities of the 
nearby bodies. By summing the eccentricities of all the plan- 
ets it is possible to average over the secular oscillations. The 
smoothed sum of the planet eccentricities shown in the top 
subpanel of Figure [4^ shows locations where stronger jumps 
in eccentricity of the entire system occur and these corre- 
spond to times when first order resonances are affecting the 
system. We can interpret the slow increases in planet eccen- 
tricity during the integration as due to first order three- 
body resonances. This follows as the zero-th order reso- 
nances should not affect the eccentricities and there are 
no strong nearby two-body resonances. As the system wan- 
ders in semi-major axis first order three-body resonances are 
crossed leading on average to the slow eccentricity evolution 
evident in Figure 2^. 



5 RESONANCE OVERLAP CRITERION 

We consider whether the number density and widths of 
three-body resonances are sufficient that they are likely to 
overlap. We estimate the number density of three body reso- 
nances and multiply this by the width of the resonances. The 
result is a filling factor that if greater than 1 implies that 
the three-body resonances overlap and so can induce chaotic 
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Figure 3. An example of a numerical integration of 5 equal mass 
tightly packed bodies experiencing three body resonance cross- 
ings. For this simulation the mass ratio m = 10~^ and initial 
interplanetary spacing 5 = 0.11247. a) The semi-major axes as 
a function of time in rotation periods of the innermost body are 
shown for all 5 bodies. Each set of points has been shifted by an 
arbitrary amount but has not been rescaled. b) We show Laplace 
angles <j) = pXi — (p + -f- qX^ in degrees as a function of 
time for p = 5,p + q = 11, q = 6 (or z = (5, —11, 6) for the in- 
ner three consecutive bodies (bottom panel; with planet indices 
i = l,j = 2,k = 3), the middle three bodies (middle panel; 
i = 2,j = 3,fc = 4) and the outer three consecutive bodies (top 
panel), c) Similar to b) except for resonances with p = 6, q = 7 (or 
z = (6, —13, 7). When one of Laplace angle ceases to circulate and 
librates about or tt, variations in the semi-major axis in three of 
the bodies are related by two conserved quantities (equation 1281 
I29I I. The middle planet moves opposite to the outer two planets 
and the middle planet experiences larger variations in semi-major 
axis than the other two. 
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Figure 4. Eccentricity evolution of the numerical integration 
shown in Figure [S] a) The top subpanel shows the smoothed sum 
of the eccentricities as a function of time. The remaining 5 sub- 
panels show the eccentricities of the 5 planets as a function of 
time, b) We show the angles iji = z ■ A in degrees as a function 
of time for z = (9, —14, 4) for the inner three consecutive bodies 
(bottom panel), the middle three bodies (middle panel), and the 
outer three consecutive bodies (top panel), c) Similar to b) except 
for z = (2, —10, 9). When the system passes through a first order 
three-body resonance variations in planet eccentricities are seen. 



behavior in the system. The 'resonance overlap criterion' 
for the onse t of chaotic behavior was pioneered by Chirikov 
in 19 59 fsee lChiriko^ll959l . [l979l : lLichtenberg fc LiebermanI 
Il992h . A similar approach has been used to estimate 
the width of t he chaotic zon e near a planet's corota- 
tion resonance ('Wisdoin| Il980l : iMurrav fc HolmanI 1 19971 : 
iQuille n & Fabcr 2006) and the onset of chaotic behavior 
in other settings (Ichirikovl [l 979: Hol man fc Murray! 
iMudrvk fc Wull2006l : lMardlinai200a ). Three-body resonance 



overlap has been seen in numerical studies of the outer Solar 
system l|Guzzdl2005l '). 

Each zero-th order three body resonance is specified by 
two integers p, q and three consecutive planets. Given a par- 
ticular p value we can first consider the separation between 
the three-body resonances with different q values. We con- 
sider an equally spaced system with separation set by 5. Let 
y be the ratio of mean motions between consecutive plan- 
ets y = {I + S)~^''^. The vector of mean motions for three 
consecutive planets n = (1, y, y^) and distance to resonance 
B = n-z (eQuations l33l[34)) . This gives B = p—{p+q)y+qy^ . 
Solving for y when B = we find that on resonance y = p/q. 
We can differentiate this with respect to q to find the dis- 
tance between resonances. Given a value for p, the distance 
between resonances with adjacent q values (in terms of dif- 
ferences in mean motions or in terms of differences in 5) is of 
order dS = pjc^ ■ For small 5, y is near 1 and so on resonance 
p ~ q. Consequently the number density of three body res- 
onances with p is (and estimated from the separation d5) 
is 



(56) 



The number density is approximately in units of the mean 
motion of the innermost body involved in the resonance. 

We now consider the width of each resonance. Consider 
5 = Spq + X where Spq is on resonance with z • n{5pq) — 0. 
We estimate the distance to resonance (equations 1331 134|l 

3 

B^n{Spq + x)-zr^ -{p~q)x (57) 



To be near or within the resonance (determined by the 
condition given in equation 21] or distance to resonance is 
smaller than y/2 times the libration frequency) 

3, 



2' 
or 



q)x\ < QmpSpq I ln<5pql exp(-(5p,p). 



(58) 



4m 



I ln<5pq| CXp( — (5p 



Here we have used equation H48p for the libration frequency. 
Multiplying by a factor of two so as to cover both sizes of 
resonance and assuming \p — q\ ~ 1 appropriate for reso- 
nances when 5 is small, the resonant width (corresponding 
to a range in 5 or 2\x\) is 



ws{p) ~ 9>mp5pq\ \n5pq \ exp(-5pgp). 



(59) 



We combine the resonant width with the number den- 
sity to estimate a three-body resonance filling factor. For 
each p the number density of resonances times their width 
is 



ps{p)ws{p) ~ Smp'^S ^llnSl exp(— (5p). 



(60) 



To estimate the total filling fraction, fs, of three-body reso- 
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nances we integrate the previous expression over all possible 
p values; 



/3 



ps{p)wsdp 



p=i 
8m5^^| In (5 



(61) 



When fs > 1 the zero-th order three-body resonances are 
sufficiently numerous and wide that they are likely to over- 
lap. 

For the simulation shown in Figure [3] we compute /a ~ 
0.75 placing the system near the regime of resonance overlap. 
This is perhaps not surprising as we found that the system 
was influenced by both the z = (6,-13,7) and (5,-11,6) 
resonances. Inverting the above equation with fa — 1 we 
find that resonance overlap occurs when 



5 < 2m 



1/4 



(62) 



The m^^'^ form of the criterion may be related to the 
slopes for stab il ity o r crossing timescales measured by 
IChambers et all (| 19961 ). 

Were we to take into account the possible combinations 
of consecutive planets (e.g, for N — 5 there are three groups 
of consecutive planets) and first order resonances, the mod- 
ified filling factor would be somewhat larger than computed 
in equation (|61|l . We expect the first order resonances to 
initially be weaker than the zero-th order resonances but 
because of the additional possible angle combinations there 
there are few times more of them. Consequently first order 
resonances may contribute to the overlap criterion. A more 
accurate criterion would likely cover the regime integrated 
by numerical stu dies where the separati on ranges from 0.5 
- 4 xm^'''' (e.g., [chambers et al.l 119961 ). The exponential 
dependence in the Laplace coefficient on interplanetary dis- 
tance (equation nop implies that three-body resonances for 
non-consecutive combinations of planets are unlikely to be 
strong. This explains why the crossing or stability timescale 
for equidistant closely spaced systems is only weakly de- 
pendent on the numbe r of planets when the number N > 5 
l|Chambers et al.1 19961 ) . The resonance overlap criterion sug- 
gests that there is a critical separation value that separates 
two regimes, an inner one at small 5 governed by instability 
from overlapping three-body resonances and an outer one 
that is much more stable. Two separate regimes and a tran- 
sition from one to another at larger sepa rations does seem 
to be exhibited in numerical integrations l|Smith fc Lissaueil 
I2OO9I ). As there are fewer combinations of consecutive plan- 
ets for low TV we would expect this transition would occur 
at smaller separatio ns; this too is seen in nu merical integra- 
tions (Figure 13 bv lSmith fc Lissauerl [20091 ). 



6 CRUDE ESTIMATES FOR DIFFUSION 

We expect that the resonances that overlap would often be 
similar in size, as illustrated in the integration shown in 
Figure [3] where the p — 5,q — 6 and p — 6,q = 7 three-body 
resonances were both important. Assuming full resonance 
overlap, we estimate a diffusion coefficient for wander in the 
Poincare coordinate 



where we use equation H38|) for the changes in A and have 
assumed that these changes take place on a timescale equal 
to the libration frequency f equation 14011. Thi s type of esti- 
mate is similar to those explored bv lChirikovl (| 19791 ). Using 
approximations for these quantities (equations |48] and I33|l 
we estimate 

Da ~ 2pm^5'^\\n5f. 

The diffusion coefficient is largest for the highest p value 
which has p ~ because the Laplace coefficients drop 
exponential at higher p. Setting p ~ 



Da ~ 2m^S-'^\lnS\'^ 



(64) 



corresponding to a diffusion coefficient in semi-major axis of 

(65) 



Da '-Sm^'^^^lhi^l^ 



An upper limit for a crossing timescale would be the 
time it takes for the semi-major axis to wander a distance 
approximate equal to the interplanetary spacing S or 

tu^S^/Da. (66) 

Using our estimate for the diffusion coefficient this gives 

(67) 



For the system we show in Figures [3] and |4] this corre- 
sponds to ~ 3 X 10^ years, and about 3 orders of mag- 
nitude gre ater than exp ected from the fit to the crossing 
times cales jFaber fc Quillcn 2007) or inferred from Figure 
3 by [chambers et aLl i 19961 ) . Consequently this timescale 
severely overestimates the crossing timescales. 

It is likely the eccentricity evolution must be considered 
as even small increases in eccentr icity can stron gly affect 
the stability or crossing timescales (|Zhou et al. \ \20m . There 
are many first order three-body resonances but eccentricity 
increases due to them are very small. There are fewer two- 
body first order resonances but eccentricity increases due 
to them could be high if the system wanders into one. The 
minimum eccentricity of a body in the vicinity o f a first 
order mean motio n resonance scales with m^^^ fe.g.. lQuillenl 
|2006| . table 1, or iMustill fc WvattI [2OI1I ) . This is a smaU 
power, and so not small for the regime covered by numerical 
integrations that have interplanetary separations of order 
a few to a dozen Hill radii. We adopt the ansatz that the 
crossing timescale is set by the time it takes for the system 
to wander into a first order mean motion resonance among 
two consecutive bodies. 

We first estimate the distance that one planet must 
wander (due to the three-body resonances) before it encoun- 
ters a first order two-body mean motion resonance. In the 
limit of high q, first order resonances (with mean motions 
with a ratio of q : g — 1) are separated in semi-major axis by 
Aa ~ q~^. For an interplanetary separation of S the nearest 
first order two-body resonance likely has q ^ . Thus the 
distance between resonances is Aa ~ . 

Using our above estimated diffusion coefficient for semi- 
major axis wander (equation[65}, the time it takes to cross a 
first order two-body mean motion resonance due to wander 
in semi- major axis would be of order 



Da ~ (AA)^a;p 



1/2 



(63) 



(68) 
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In Figure [5] we show a comparison between crossing times 
predicted with the above t2 timescale and crossing times 
measured numerically. The numerically measured timescales 
are shown with the fit t o the numerically measured value by 
iFaber fc QuiiienI l|2007l ). This estimate is within two orders 
of magnitude of the relation found numerically and overesti- 
mates the crossing time at small separations. Nevertheless it 
is the first analytical derived estimate of crossing time that 
covers an appropriate range in parameter space. We note 
that because of the high power of S in the above equation 
our power law relation is nearly as steep as the numerically 
measured times that have primarily been fit with exponen- 
tial functions. On Figure [S] the resonance overlap boundary 
would be at a constant value of S/jj}'''^ corresponding to a 
vertical line on the right hand side of this plot. Numerical 
integrations cover the range of S/fi^^"^ shown in this plot 
hence we expect that the overlap criterion line should lie 
on the right hand side of the plot with 5/^^^* ~ 3.5. Our 
predicted location (equation I62|l for the onset of three-body 
resonance overlap is about a factor of two or so too low as 
this line lies in the middle rather than the right hand side 
of the plot. 

The above crude estimates for a diffusion coefficient 
(equation I64|) and a crossing time (equation [68} are strong 
power law functions of both mass and interplanetary sep- 
aration. Their strong dependence on separation suggests 
that the exponential forms fit to the num erically mea- 



sured stability or crossing timescales (Chambc rs~et al.l 19961: 
Duncan & Lissauer 1997: 'Z hou et al.ir2007l : lChatteriee et al.l 
2008; Smith fc Lissauer ,20091 )" might in future be accounted 
for through chaotic motions induced by three-body reso- 
nances. Perhaps the exponential form has provided a good 
fit because of the limited range of times over which these 
systems can be integrated and because the diffusion rate is 
such a strong function of mass and separation. If three-body 
resonances overlap then there is no underlying mathemati- 
cal reason (related to Arnold diffusion or the Nekhoroshev 
theorem) that would predict an exponential dependence on 
interplanetary separation and mass. We suspect that only 
outside the regime of three-body resonance overlap could a 
true long timescale exponential dependence on planet mass 
and separation be recovered. 



7 SUMMARY AND DISCUSSION 

In this paper we have considered the role of three-body 
resonances for an idealized equal mass, uniformly spaced, 
but closely packed initially low eccentricity co-planar mul- 
tiple planet system. We have estimated the strengths and 
libration timescales for zero-th order (in eccentricity) three- 
body resonances using an asymptotic approximation to the 
s = 1/2 Laplace coefficient. Two conserved quantities are 
present relating variations in semi-major axes between the 
bodies affected by the resonance. These variations and the 
Laplace angles are useful for identifying the effect of three- 
body resonances in numerical integrations of multiple planet 
systems. 

By estimating the number and widths of the three- 
body resonances, we have derived an approximate reso- 
nance overlap criterion. We find that zero-th order three- 
body resonances are likely to overlap when the separa- 
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Figure 5. A comparison of numerical to estimated crossing 
timescales. The crossing time estimated from the time it takes 
the system to cross a first order mean motion resonance between 
two consecutive bodies (calculated with equation I68p is shown 
as points for three different planet mass ratios. The function 
fit to numericall y mea sured crossing times (using equation 2 by 
iFabcr fc Quilleni l2007h is shown as line segments for the same 
three planet mass ratios. 



tion between planets S < 2m^''*. The resonance over- 
lap criterion is close to the regime covered by numeri- 
cal integr ations of multiple planet systems exhibiting in- 
stability (IChambers et ahl 19961; Duncan fc Lissauer! 19971; 
Zhou et al] l2007l : IChatteriee et al.l ,2008, : , Smith fc Lissaueij 



20091 ). suggesting that the instability seen in integrations 
of closely spaced multiple planetary systems is due to 
chaotic behavior associated with multiple three-body res- 
onances. We note that previous studies of asteroids have 
also attributed chaotic behavior to three-body r esonances 
(|Murrav et al.lll998l : iNesvornv fc Morbidellilll998D . 

Our resonance overlap criterion lies near but within the 
regime covered by numerical integrations exhibiting instabil- 
ity implying that we have underestimated the filling factor of 
resonances by a factor of a few. However, we have not taken 
into account first order three-body resonances, the different 
combinations of consecutive planets, indirect terms in the 
Hamiltonian and we have only crudely estimated resonance 
strengths. Future works can improve upon the overlap cri- 
terion by expanding and improving the calculation. 

For spacings larger than the overlap criterion boundary 
three-body resonances should not be as dense and the prob- 
ability of resonance overlap drops. We postulate that there 
is a region of long timescale stability at large separations. 
This region and the transition between the two regimes is 
likely the reason for measured changes in slope of crossing 
time versus separation and a strong increa se in crossing time 
at la rge separations (see Figures 1-4 by ISmith fc Lissauer! 
!2009D . The filhng factor of three-body resonances should be 
lower for systems with fewer planets because there are fewer 
combinations of consecutive planets and so fewer strong 
three-body resonances. Thus we expect the transition to a 
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more stable regime would occur at smaller planetary sep- 
arations when there are fewer planets. This also is seen 
in numerical integratio ns (see red points in Figure 13 by 
ISmith fc LissaueIj|2009^ . 

We have attempted to predict diffusion rates using 
three-body resonances. The timescale to wander a distance 
of the interplanetary separation grossly overestimates the 
crossing timescale whereas that to diffuse until the system 
crosses a first order mean motion resonance between two 
bodies overestimates the crossing timescale by 1 or 2 orders 
of magnitude at small separations. This estimated timescale 
depends on the 8-th power of the interplanetary spacing 
suggesting that exponential functions have primarily been 
successful at fitting numerically measured crossing timescale 
because of the strong dependence on separation of the three- 
body resonances. Future work could strive to improve upon 
these gross estimates. Crossing timescales are not directly 
related to diffusion coefficients and the dy namics may be in- 
termittent a nd diffusion anisotropic (e.g.. IShevchenkoll2010l : 
lGuzzoll2005l ). To better account for or predict the crossing 
timescales perhaps Lyapunov timescales and diffusion coef- 
ficients (i.e., eccentricity growth rates and rates of wander in 
semi-major axis) could be measured directly from numerical 
integrations. These then may be easier to understand with 
analytical estimates such as explored here. 

Here we have considered equal mass, equidistant copla- 
nar systems. However much of the framework developed 
here could be applied to less ideal systems such as closely 
spaced multiple planet extrasolar planetary systems. We re- 
mind the reader that here we have focused on systems that 
are in three-body resonances but are not in strong two- 
body resonances. Three body resonances may also be im- 
portant in these systems but calculations are likely to be 
more challenging in this setting. Diffusion in semi-major 
axis seen in simulations of closely spaced satellite systems , 
e.g., the Uranian satellite system. (iDuncan fc Lissaueill 19971 : 
IShowalter fc Lissau ei"200ffi'Dawson et al.lbOld ') and the Ke- 
pler 11 system (|Lissauor ct al. 2010) might in future be in- 
terpreted in terms of variations arising from three-body res- 
onances. 
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